Resurrecting extinct cephalopods with biomimetic robots to explore hydrodynamic stability, maneuverability, and physical constraints on life habits

Externally shelled cephalopods with coiled, planispiral conchs were ecologically successful for hundreds of millions of years. These animals displayed remarkable morphological disparity, reflecting comparable differences in physical properties that would have constrained their life habits and ecological roles. To investigate these constraints, self-propelling, neutrally buoyant, biomimetic robots were 3D-printed for four disparate morphologies. These robots were engineered to assume orientations computed from virtual hydrostatic simulations while producing Nautilus-like thrusts. Compressed morphotypes had improved hydrodynamic stability (coasting efficiency) and experienced lower drag while jetting backwards. However, inflated morphotypes had improved maneuverability while rotating about the vertical axis. These differences highlight an inescapable physical tradeoff between hydrodynamic stability and yaw maneuverability, illuminating different functional advantages and life-habit constraints across the cephalopod morphospace. This tradeoff reveals there is no single optimum conch morphology, and elucidates the success and iterative evolution of disparate morphologies through deep time, including non-streamlined forms.


Supplementary Text
Measurements of septal spacing were recorded from a CT-scanned Nautilus pompilius conch. These measurements were recorded as the angle between the ventral edge of the current and previous septum and the spiraling axis of the shell. Because septal spacing differs in early ontogeny (Fig. S11), only measurements from the 7 th to 33 rd (terminal) septum were considered. The average angle of 23.46° ± 3.32° (standard deviation) was rounded to 23° and held constant through the ontogeny of the hydrostatic models.
Measurements of shell and septum thickness were recorded as a ratio of inner whorl height (Table S13) in order to apply thickness to the theoretical planispiral models. These models were created with array instructions (Table S14) in Blender (Blender Online Community) to build the shell from measurements of representative specimens (Sphenodiscusoxycone; Dactyliocerasserpenticone; and Goniatitessphaerocone).
The whorl section of the center model is an average (Fig. S12) of the near-endmember models (oxycone, serpenticone, and sphaerocone).
Westermann morphospace parameters were computed from measurements on each virtual model ( Fig. S13; Table S15). These parameters fall on a ternary diagram with variable whorl expansion (w), umbilical exposure (U), and inflation (i.e., thickness ratio; Th). Each are computed with the equation below with variables corresponding to those depicted in Figure S13: w = a/a' (S1) U = UD/D (S2) Th = b/D (S3) The peripheral shape of the Nautilus septum (suture line) and its general first-order morphology was preserved while modifying it to fit within the whorl sections of the theoretical models (Fig. S14).

Calculation of drag coefficients (Cd)
Drag coefficients do not capture the hydrodynamics of a single shape moving at different velocities (or with changing scale; i.e., Reynolds numbers). However, these coefficients can be used to approximate the relative differences in hydrodynamic drag between different morphologies. Drag coefficients were computed from analyzing the deceleration of each biomimetic robot during the one-second pulse experiments after they ceased jetting (i.e., as they fell from their maximum attained velocities). These coefficients were determined in MATLAB by creating a modeled velocity function that reduced the sum of squared error between this function and the data points computed from all trials for each examined morphology. The velocity function (Equation S5) is derived from the drag force equation (Equation S4) and follows the form: Where u(t) is the velocity as a function of time (t), m is mass, u0 is the initial velocity, V is volume (note that this becomes equivalent to the area term, and that volume is nearly identical between each robot), ρ is density, Cd is the drag coefficient. Drag coefficients were fit with the curve fitting toolbox in MATLAB. Note that the time term in the equation has 1 second subtracted. This value corrects for the time the motor was active and translates the curve to the moment the robots stopped jetting. Drag coefficients and modeled velocity function terms are reported in Fig. S6 and Table S5.

Calculation of moments of inertia
Moments of inertia for each material of unique density for the robots (PETG thermoplastic, electronics cartridge, bismuth counterweight, pump and chamber liquid, motor (bulk ρ), both batteries (bulk ρ), electronics (bulk ρ), and self-healing rubber) and the virtual hydrostatic models representing the living animals (soft body and shell) were computed to determine their total moments of inertia (Table S10; sum of each component). The moments of inertia for each individual component were computed with the program MeshLab. The relative importance of the moments of inertia versus hydrodynamic effects (drag and wake dynamics, etc.) were determined by comparing the angular velocities computed for rotation in a vacuum (after one second resulting from 0.3N of thrust) and those observed after one second of jetting in the robots. Angular velocities (ω) in a vacuum were computed with the following equations: Where is torque, F is the force of thrust (~0.3N for the robots), r is the length of the lever arm measured between the source of jet thrust (i.e., the hyponome) and the vertical axis passing through the centers of buoyancy and mass (measured from the virtual models in the program Blender; Table S11), I is the total moment of inertia (Table S10), and α is angular acceleration. ωf and ωi are angular velocities after one second of jetting and just before jetting (equal to zero), respectively, and Δt is the change in time (equal to one second). Potential differences in rotational kinematics due to different moments of inertia between the robots and the virtual hydrostatic models (representing theoretical morphologies of the living animals) were estimated by comparing the proportion of the observed ω values after one second of jetting, and the ω values computed from each total moment of inertia value (Table S11). Figure S1: Schematic of the thrust calibration setup. Each robot was placed in water with a fishing line attached from the hyponome location, through a series of pulleys, to a dual-range force sensor. Rendered in MeshLab.

Figure S2
: Thrust calibration of biomimetic cephalopod robots. A) Example of a single recorded trial. Force was measured for 30 seconds to capture the 6-second region of motor activity (highlighted grey). The model pulled a fishing line on a pulley system, attached to the force sensor. When the line became taught, force peaked, then oscillated. Only the stable region (denoted between the two dashed lines) was used to measure the average thrust for each trial. After isolating this region, the true zero-datum was subtracted from the data series. B) Average thrust recorded during 15 trials for each model with the motor running at 100% maximum voltage (7.4V). C) Average thrust recorded during 15 trials for each model with voltage adjusted via pulse-width modulation (Serpenticone = 100%; Oxycone = 100%; Sphaerocone = 95%; Center = 85%. Error bars represent one standard deviation. Figure S3: Rocking during movement demonstrated by computing the angle displaced from the static orientation. A) Serpenticone, B) oxycone, C) sphaerocone, and D) morphospace center. These models with artificially high hydrostatic stability, experience minimal rocking (+/-5 degrees). Figure S4: Example frame of the motion tracking footage showing grabber tool used to position the biomimetic robot. After releasing the robot, a fiber optic cable was used to deliver an infrared pulse that initated the motor. At the same time, a green indicator LED illuminated the model while the motor was active. This light was used to determine time-zero for each experiment. Various kinematics were monitored with the two tracking points placed on each model. Figure S5: Example of drag coefficient and Reynolds number curves for conch shapes similar to the cephalopod robots. Since drag coefficients vary based on size and swimming speed, hydrodynamic properties are interpreted in terms of velocity and acceleration. Figure S6: Analysis of drag coefficients for the serpenticone (A), oxycone (B), sphaerocone (C), and morphospace center (D). Drag coefficients (E) were determined by minimizing the sum of squared errors between the modeled velocity function (black) and all computed datapoints for each morphotype. The errors bars on panel E denote 95% confidence intervals for the drag coefficients. Note that the serpenticone and sphaerocone cannot be statistically distinguished.   Figure S10: Average duration of all trials for the serpenticone (n = 17), oxycone (n = 10), sphaerocone (n = 12), and morphospace center (n = 9). Error bars represent standard deviations. A one-way ANOVA with a Games-Howell post hoc test was used to distinguish means. Note that the differences of all other combinations than those denoted are significant at the p ≤ 0.0001 level.
Figure S11: Septal spacing for each septum recorded from a CT-scanned conch of Nautilus pompilius. The shaded region corresponds to the angles that were considered for the theoretical conch models (7 th to 33rd septum).        Table S5: Drag coefficients (Cd) computed for each robot by analyzing the deceleration of all corresponding trials after they ceased jetting (i.e., falling from their maximum attained velocities). Each term in the modeled velocity function is reported below (m = mass, u0 = initial velocity the moment jetting ceased, ρ = water density, V = volume which is raised to the 2/3 power to equal the area term, Cd = drag coefficient, Adj. r 2 = adjusted r-squared. Note that 95% upper and lower confidence intervals (CI+ and CI-, respectively) are listed for each computed drag coefficient. The serpenticone and center robots have overlapping confidence intervals, suggesting they cannot be statistically distinguished from each other.      Table S10: Moments of inertia (MOI) computed for each component of unique density (including some bulk density values; see methods) for the robots, virtual models representing living animals with theoretical morphologies, and two simple shapes ("Sphere" and "Disk"). The sphere has the same volume as the sphaerocone, and the disk has the same aspect ratio as the oxycone. All MOIs were computed in MeshLab. The simple shapes were calculated by hand for comparison.   (Table S10) for the robots and hydrostatic models (representing the living animals) were used to compute angular velocity ω reached after one second of jetting at 0.3N (rows b and c, respectively). The observed ω values after jetting for one second (row d) are consistently about an order of magnitude lower, demonstrating that hydrodynamic effects dominate rotational kinematics. The hydrodynamic proportions of ω for the robot moment of inertia and "living animal" moment of inertia (rows e and f, respectively) were subtracted to determine the relative difference in these proportions (row g). This approach yields estimates of how much the lower moments of inertia in the robots influence rotational kinematics for each morphotype.       Vct (cm 3 ) msb (g) msh (g) mcl (g) mcg (g) mwd (g) mtotal (g)